
function [z, res] = ALSbaseline(y, lambda, p, itter)
% Estimate baseline with asymmetric least squares
% 
% Reference:
% Eilers, P. and H. Boelens "Baseline Correction with Asymmetric Least
% Squares Smoothing."
% http://www.science.uva.nl/~hboelens/publications/draftpub/
%
% INPUT:    y       = Signal
%           lambda  = smoothness     
%                           "and 10^2 > lambda > 10^9, but exceptions
%                               may occur." 
%           p       = asymmetry      
%                           "We found that generally 0.001 > p > 0.1 is a
%                               good choice (for a signal with positive
%                               peaks)
%           itter   = number of interratiosn    (DEFAULT = 10)
%
% OUTPUT:   z       = calculated baseline
%           res     = residuals (deviation of baseline from signal)

% 

if nargin<4,
    help ALSbaseline
    return
end

if itter < 1 || itter > 1000
    ii = 10;
else ii = itter;
end
    
m = length(y);
D = diff(speye(m), 2);
w = ones(m, 1);

for it = 1:ii
W = spdiags(w, 0, m, m);
C = chol(W + lambda * D' * D);
z = C \ (C' \ (w .* y));
w = p * (y > z) + (1 - p) * (y < z);
end
    res = y - z;
